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FIELD OF THE INVENTION 

The present invention relates to a method for simulating well tests on a discrete 
fractured reservoir model. 

The method can be implemented for example in the field of oil production by 
5 reservoir engineers in order to obtain reliable flow predictions. 

A well test is to be simulated in a permeable porous medium (reservoir) crossed 
through by a network of fractures much more conducting than the porous matrix. 
Fractured reservoirs constitute an extreme type of heterogeneous reservoirs comprising 
two very different media, a matrix medium containing the most part of the oil in place 
10 and having a low permeability, and a fractured medium representing less than 1 % of 
the oil in place and highly conducting. The fractured medium itself can be complex, 
with different series of fractures characterized by their respective density, length, 
orientation, inclination and opening. 

During well testing, the flow rate conditions imposed on the well lead the oil 
15 contained in the reservoir to flow towards the well. It is a single-phase (the only mobile 
phase is the oil) and compressible flow (the rock of the reservoir and the oil fluid are 
compressible). 

For any elementary volume of the reservoir, the pressure of the oil contained in this 
volume is governed by the following equation, if gravity is not taken into account : 



(P.C T .^- = div(— NP) + Q (1) 
d t /J 
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where : 

(|) : represents the pore volume, 
C T , the total compressibility (fluid + rock), 
K, the permeability of the rock, 
5 u, the viscosity of the fluid, 

Q, the incoming flow, which is zero everywhere except in places where the well 
communicates with the reservoir, and 
P, the unknown pressure. 

In order to simulate a well test, whatever the medium, this equation has to be solved 
10 in space and in time. Discretization of the reservoir (mesh pattern) is therefore 
performed and solution of the problem consists in finding the pressures of the meshes 
with time, itself discretized in a certain number of time intervals. 

In the case of a fractured medium, if all the complexity of the fracture network is to 
be taken into account, it is necessary to have a mesh pattern that represents this network 
15 accurately. Furthermore, in such a medium, the permeabilities K are generally much 
higher in the fractures than in the matrix, and the flows are consequently fast in the 
fractures and slow in the matrix. 

BACKGROUND OF THE INVENTION 

Single-phase flow computing tools are currently used; however, they are not applied 
20 to the real geologic medium in all its complexity but to a homogenized representation, 
according to the reservoir model called double-medium model described for example by 
Warren and Root in "The Behavior of Naturally Fractured Reservoirs", SPE Journal, 
September 1963. Any elementary volume of the fractured reservoir is thus modelled in 
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the form of a set of identical parallelepipedic blocks limited by an orthogonal system of 
continuous uniform fractures oriented in the direction of one of the three main 
directions of flow (model referred to as "sugar box" model). The fluid flow on the 
reservoir scale occurs through the fractured medium only, and fluid exchanges take 
5 place locally between the fractures and the matrix blocks are not applied to the complex 
medium. This representation, which does not reproduce the complexity of the fracture 
network in a reservoir, is however effective but at the level of a reservoir mesh whose 
typical dimensions are 100 m x 100 m. 

It is however preferable to carry out the flow calculations on the "real" geologic 
10 model instead of an equivalent homogeneous model, which affords the double 
advantage of allowing : 

- to validate the geologic image of the reservoir built by the geologist from all the 
information he has gathered on the reservoir fracturation (this validation being 
performed by comparison with the well test data), 

15 - to reliably test the sensitivity of the hydraulic behaviour of the medium to the 
uncertainties on the geologic image of the fractured medium. 

A well-known modelling method consists in finely meshing the fracture network 
and the matrix while making no approximation concerning fluid exchanges between the 
two media. It is however difficult to implement because the often complex geometry of 
20 the spaces between the fractures makes it difficult to mesh them and, in any case, the 
number of meshes to be processed is finally often tremendous. The complexity 
increases still further with a 3D mesh pattern. 
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Techniques for modelling fractured porous media are described in patents FR-A- 
2,757,947 and FR-A-2,757,957 filed by the applicant. The former technique relates to 
determination of the equivalent fracture permeability of a fracture network in an 
underground multilayer medium from a known representation of this network, allowing 
5 to systematically connect fractured reservoir characterization models to double-porosity 
simulators in order to obtain more realistic modelling of a fractured underground 
geologic structure. 

The second technique relates to simplified modelling of a porous heterogeneous 
geologic medium (such as a reservoir crossed through by an irregular network of 
10 fractures for example) in the form of a transposed or equivalent medium so that the 
transposed medium is equivalent to the original medium, according to a determined type 
of physical transfer function (known for the transposed medium). 

SUMMARY OF THE INVENTION 

The method according to the invention allows to simulate single-phase flows in a 
15 fractured medium. 

The object of the method according to the invention is to model fluid flows in a 
fractured multilayer porous medium by taking account of the real geometry of the 
fracture network and the local exchanges between the porous matrix and the fractures at 
each node of the network, and thus to allow simulation of the interactions between the 
20 pressure and flow rate variations in a well running across said medium. It is 
characterized in that the fractured medium is discretized by means of a mesh pattern 
wherein the fracture meshes are centred on nodes at the various intersections of the 
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fractures, each node being associated with a matrix volume, and the flows between each 
fracture mesh and the associated matrix volume are determined in a pseudosteady state. 

The matrix volume associated with each fracture mesh in each layer is for example 
delimited by all of the points that are closer to the corresponding node than to 
5 neighbouring nodes. 

In order to determine the matrix volume associated with each fracture mesh, each 
fractured layer is for example discretized in pixels and the distance from each pixel to 
the closest fracture mesh is determined. 

The value of the transmissivity for each fracture mesh - matrix block pair is 
10 determined for example by considering that the pressure varies linearly as a function of 
the distance from the point considered to the fracture mesh associated with the block. 

The method according to the invention greatly simplifies calculation of the 
exchanges between the matrix and the fractures. In fact, there is no specific mesh 
pattern for the matrix, but each node of the fracture network is associated with a matrix 
15 volume, and the flow between each fracture node and its associated block is calculated 
in a pseudosteady state. The innovation concerns calculation of the volume of each 
block and of the proportionality coefficient (called transmissivity) that connects the 
matrix-fracture flow rate to the matrix-fracture pressure difference. 

The main advantage of the method according to the invention in relation to methods 
20 applied to a finely meshed real fracture network lies, as detailed in the description 
hereafter, in the considerable reduction in the number of meshes centred on fracture 
nodes used to solve equation (1). It can be checked that the simulation rate saving in 
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relation to prior fine-mesh methods is greater than the simple arithmetical ratio of the 
number of meshes used in both cases. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Other features and advantages of the method according to the invention will be 
5 clear from reading the description hereafter of a non limitative realisation example, with 
reference to the accompanying drawings wherein : 

- Figure 1 shows an example of a 2D fracture mesh, FM, 

- Figure 2 shows an example of a 2D matrix block, MB, 

- Figure 3 shows an imposed flow rate variation curve F(t) in a well test, and 

10 - Figure 4 shows an example of a fracture network for which the method according to 
the invention leads to a radical reduction in the number of meshes to be processed. 

DETAILED DESCRIPTION OF THE METHOD 

1) Fracture network discretization 

The method comprises meshing the multilayer fractured reservoir modelled for 
15 example by means of the fractured reservoir modelling technique described notably in 
the aforementioned patent FR-2,757,947, assuming that the fractures are substantially 
perpendicular to the layer boundaries. 

In order to form this mesh pattern, all the fracture intersections or nodes CN are 
located (Fig.l) and fracture meshes are formed by associating a single matrix block with 
20 each node. For the purpose of 3D modelling, additional nodes have to be added in the 
neighbouring layers in order to take account of the fractures that run across several of 
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them. The computing nodes thus defined constitute the centre of the fracture meshes 
FM. The meshes are given boundaries such as the ends of the fractures on the one hand 
and the midpoint of the segments connecting the computing nodes on the other hand. 
Considering these boundaries imposed on the meshes, their volume <j) can be calculated 
5 (see relation 1). Calculation of the connections between fracture meshes 
(transmissivities) used for calculation of the inter-mesh flows can be carried out 
according to the method described in the aforementioned patent FR-2,757,947. 

2) Matrix medium discretization 

Discretization of the matrix medium consists in assigning a rock volume to each 
10 fracture mesh. During dynamic simulation, exchanges between the matrix and the 
fractures will be calculated in the pseudosteady state (exchange flow proportional to the 
pressure difference) between each fracture mesh and its associated single matrix block. 

For calculation of these matrix volumes, the problem is dealt with layer by layer, 
which amounts to disregarding the vertical flows in relation to the horizontal flows in 
15 the matrix. 

For a given layer, the blocks are defined as follows : 

- the block heights are equal and fixed to the height of the layer, 

- in the plane of the layer, the surface area of the block associated with a fracture 
mesh is all the points that are closer to this fracture mesh than to another one. 

20 Physically, this definition implies that the boundaries at zero flow in the matrix are 

the equidistance lines between the fracture meshes. This approximation is valid if the 
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neighbouring fractures are at very close pressures, which is true in the presence of 
highly conducting fractures in relation to the matrix. 

a) Block geometry 

In order to determine the geometry of the blocks in a given layer, a two-dimensional 
5 problem consisting in finding, for each fracture mesh, the points that are closer to this 
mesh than to another one has to be solved. This problem is solved by using for example, 
but preferably, the geometric method described in the other aforementioned patent FR- 
2,757,957, which allows, by discretizing the fractured layer into a set of pixels and by 
applying an image processing algorithm, to determine the distance from each pixel to 
10 the closest fracture. During the initialization phase of this algorithm, value 0 (zero 
distance) is assigned to the pixels belonging to a fracture and a high value is assigned to 
the others. Furthermore, if the number of the fracture mesh to which each fracture pixel 
belongs is given in this phase, the same algorithm finally allows to determine, for each 
pixel : 

15 - the distance from this pixel to the closest fracture mesh, 
- the number of the closest fracture mesh. 

All of the pixels having the same fracture mesh number constitute the surface area 
of the matrix block associated with this mesh. Multiplying this surface area by the 
height of the layer allows to obtain the volume of the block associated with the mesh. 
20 Figure 2 shows an example of a thus obtained 2D matrix block. 

By construction, the sum of the surface areas of the blocks is equal to the total 
surface area of the layer, which guarantees that no volume is added to or taken from the 
layer. 
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For each matrix block, the image processing algorithm also gives the distance from 
each pixel of the block to the associated fracture mesh. This data is used to calculate the 
transmissivity between the fracture mesh and the matrix block. 

b) Matrix-fracture transmissivity 

5 On the assumption of a pseudosteady state where the exchange flow is considered 

to be proportional to the difference in the pressures of each one of them, we have the 
following relation : 

V- 

where : 

F mf is the matrix-fracture flow, 
10 T^, the matrix-fracture transmissivity, 
u, the viscosity, 

P m , the pressure of the matrix block, and 

P f , the pressure of the associated fracture mesh. 

On the assumption of a pseudosteady state, the value of transmissivity T mf is 
15 constant during simulation. This value is here calculated for each fracture mesh - matrix 
block pair. 

For this calculation, it is assumed that, in the matrix block, the pressure varies 
linearly as a function of the distance from the point considered to the fracture mesh 
associated with the block. The transmissivity is defined as follows : 
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where : 

1 is the length of the fractures in the fracture mesh, 
H, the height of the layer (and of the matrix block), 
K, the permeability of the matrix, and 

D, the distance from the fractures for which the pressure is the mean pressure of the 
block. 

1, H and K are known (factor 2 comes from the fact that the fracture has two faces). 
Calculation of D is made from the information provided by the image processing 
algorithm concerning the distances from the pixels to the closest fractures. In fact, 
according to the hypothesis of linear variation of the pressure as a function of the 
distance to the fracture, we have : 




where S is, in 2D, the surface area of the matrix block and d(s) is the function giving the 
distances between the points of the matrix block and the associated fracture mesh. 

In terms of pixels, we thus simply have : 

where N is the number of pixels of the matrix block and d n the distance from pixel n to 
the fracture mesh. 
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3) Well representation 

A well is represented by its geometry on the one hand and by the flow rate 
constraints imposed thereon with time on the other hand. 

a) Geometry 

A well consists of a series of connected segments that intersect the network 
fractures. The geometric representation of a well is therefore a 3D broken line. Some 
segments of this broken line communicate with the reservoir. The points of 
communication between the well and the network are the points of intersection of these 
segments communicating with the fractures. 

b) Flow rate constraints 

The flow rate constraints are imposed at the point where the well enters the 
fractured medium. They are entered by the user in the form of a step function giving the 
flow rate as a function of time. The flow rate change times of the well indicate the 
various periods to be simulated. 

For a conventional well test (Fig.3), a flow rate F(t) is imposed for a first period 
after which the well is closed (zero flow rate) and the pressure buildup is observed in 
the well. 

4) Dynamic simulation 

During dynamic simulation, the simulated period of time is split up into time 
intervals whose duration ranges, for a well test, between 1 s (just after opening oi 
closing of the well) and some hours. 
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Passage from the time n (where all the pressure values are known) to the time n+l 
is achieved by solving, in all the meshes and blocks of the reservoir, the following 
equation which corresponds to relation 1 once discretized : 

^.c r . ^"T +£--^" +l - p ; +1 ) = 2 (5) 

t n+l -t n Y m 

where : 

i, the number of the mesh or of the block, 

j, the number of the mesh or of the block next to i, 

P n , the pressure at the time n, 

t n , the time at the time n, 

Q b the flow entering i, and 

Ty, the connection transmissivity between i and j. 

Apart from the pressures, the terms of this equation are known. The pore volumes 
of the fracture meshes and of the matrix blocks ($) are known by means of the mesh 
pattern, the fracture-fracture and matrix-fracture transmissivities (T,j) are calculated as 
described above and flow rates (Qi) are zero everywhere except at the well-reservoir 
connections where they are imposed. 

a) Elimination of the unknowns relative to the matrix 

Equation (6) can be solved by inversion of a linear system where the unknowns are 
the pressures of the fracture meshes and the pressures of the matrix blocks at the time 
n+l. 
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However, since each matrix mesh is only connected to one matrix block and since 
the blocks are not connected to each other, it is possible to eliminate the unknowns of 
the blocks of the system whose size is then divided by 2. 

In fact, for a matrix block, equation (6) is expressed as follows : 

5 subscripts m and f denoting the matrix and the fracture respectively, and C Tm being the 
total compressibility of the matrix. 

We have C Tm = c m + Sw.c w + (1-Sw).c 0 , where Sw is the residual water saturation in 
the matrix, c m the compressibility of the matrix, c w the compressibility of the water and 
c 0 the compressibility of the oil, and <\> m is the pore volume of the matrix block, i.e. the 
10 volume of the block multiplied by the porosity of the matrix. 

We deduce therefrom that : 



with : 

<P m -C Tm 



t n+l_ t n 



By transferring this expression into the equation relative to the fracture mesh, we 
obtain : 
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A + 



.r l -S--C = 




■Pi 



(9) 



where i is the number of the fracture mesh considered, j the numbers of the fracture 
meshes connected to i, and m the number of the block associated with fracture mesh i. 

Furthermore, we have : 



where C Tf is the total fracture compressibility : C Tf = c f + c 0 . 

Reduction of the number of unknowns by a factor 2 is a further advantage of the 
present method. The solution time for a linear system develops approximately like the 
square of the number of unknowns, and consequently elimination of the matrix 
unknowns leads to a considerable saving in the rate of simulation of a well test on a 
fractured network. 

b) Time interval solution algorithm 

The general algorithm for solving a time interval consists in forming a linear system 
corresponding to equation (9) relative to any fracture mesh i, in solving this linear 
system and in updating the unknowns. 

The stages of this algorithm are as follows : 

- Forming the linear system with calculation of the terms of equation (9) 



Calculation of the fracture medium contribution 



accumulation term (A t ) 
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connection between fracture meshes (Tj/u) 
well boundary conditions (Q;) 

Calculation of the matrix medium contribution 
accumulation term (AJ 
matrix-fracture connection term (U m ) 

- Solving the linear system (so as to find the fracture pressure values at the time n+1) 
by means of the iterative conjugate gradient algorithm (CGS) 

- Updating with calculation of the matrix block pressures using equation (8), and time 
interval management. 

c) Time interval management 

During well test simulation, the time intervals are generally very short after opening 
or closing of the well, because the pressure variations are high at these times. Later, the 
time intervals can be longer when the pressure variations are lower. 

Management of the time intervals consists in connecting the duration of the time 
interval n+1 to the maximum pressure variation observed during time interval n. The 
user must therefore provide the following data : 

Initial time interval : dt mi 
Minimum time interval : dt mn 
Maximum time interval : dt max 
Lower limit of pressure variation : dp^ 
Upper limit of pressure variation : dp max 
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Time interval increase factor : rdt + (>1) 
Time interval reduction factor : rdt - (<1). 

From these values, management of the time interval is performed as follows : 

At the time t n , the majorant of the pressure variations in the fracture meshes and the 
matrix blocks during time interval n is calculated (i.e. between times t"" 1 and t n ). 
According to the value of this pressure variation dp, the time interval is either : 

increased by a factor rdt+ if dp < dp^, 

decreased by a factor rdt- if dp > dp max , 

unchanged in the other cases. 

After each well flow rate condition change, the time interval is re-initialized to the 
value dt mi . Furthermore, an optimization is performed at the end of the simulation 
period : if tf is the time to be simulated to reach the end of the period and dt the planned 
time interval, the time interval selected is min(t f/2 ,dt). 

5) Simulation input data 

a) Fracture network 

The geometry of the fracture network and the attributes of the fractures 
(conductivity, opening) are given in the form of a file as in the method described in the 
aforementioned patent FR-2,7 57,947. 

b) Petrophysics 

The petrophysical data to be provided are as follows (the unit is given between 
parentheses) : 
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- compressibility of the fracture network (c f in bar 1 ) 

- compressibility of the matrix (c m in bar" 1 ) 
matrix permeability (K in mD) 

matrix porosity ((j) m dimensionless). 

5 The data in question are considered to be homogeneous on the fractured medium 

considered. 

c) Fluids 

The data relative to the fluids contained in the fractured medium concern the oil 
(mobile) and the water (immobile) that is always present in residual amounts in the 
10 matrix : 

- residual water saturation in the matrix (Sw dimensionless) 

- compressibility of the water (c w in bar 1 ) 
compressibility of the oil (c 0 in bar" 1 ) 

- viscosity of the oil (u in cpo). 

15 d) Well 

The data relative to the well are : 

- its geometry, in the form of a series of connected segments 

- the imposed flow rates, in the form of a curve giving the imposed flow rate as a 
function of time. 
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e) Time management 

The values to be provided are the values already defined in the chapter relative to 
the time interval management : dt mi , dt^ dt max , dp min , dp™*, rdt+ and rdt-. 

Comparative meshing example 

The following comparative example illustrates the meshing simplification and the 
correlative saving in flow calculating time that can be obtained by implementing the 
method. By way of comparison, it has been observed that the fractured zone shown in 
Fig.4, that normally would have been meshed with 80,000 meshes with conventional 
meshing techniques, can be meshed here with about 4000 meshes only. 



19 



CLAIMS 

1) A method for modelling fluid flows in a fractured multilayer porous medium in 
order to simulate interactions between pressure and flow rate variations in a well 
running across said medium, characterized in that the fractured medium is discretized 

5 by means of a mesh pattern wherein the fracture meshes are centred on nodes at the 
various fracture intersections, each node being associated with a matrix volume, and 
flows are determined between each fracture mesh and the associated matrix volume in a 
pseudosteady state. 

2) A method as claimed in claim 1, characterized in that the matrix volume 
10 associated with each fracture mesh is delimited in each layer by all of the points that are 

closer to the corresponding node than to neighbouring nodes. 

3) A method as claimed in claim 2, characterized in that each fractured layer is 
discretized in pixels and the matrix volume associated with each fracture mesh is 
delimited by determining the distance from each pixel to the closest fracture mesh. 

15 4) A method as claimed in any one of the previous claims, characterized in that the 

transmissivity value is determined for each fracture mesh - matrix block pair by 
considering that the pressure varies linearly as a function of the distance from the point 
considered to the fracture mesh associated with the block. 



PATENT 



INSTITUT FRANQAIS DU PETROLE 

METHOD FOR MODELLING FLUID FLOWS 
IN A FRACTURED MULTILAYER POROUS MEDIUM 
AND CORRELATIVE INTERACTIONS IN A PRODUCTION WELL 

Inventor : Sylvain SARDA 
ABSTRACT 

- Method for modelling fluid flows in a fractured multilayer porous medium by taking 
account of the real geometry of the fracture network and the local exchanges between 
the porous matrix and the fractures at each node of the network, thus allowing 
simulation of the interactions between the pressure and flow rate variations in a well 
running across said medium. The method essentially comprises discretization of the 
fractured medium by means of a mesh pattern, with fracture meshes centred on nodes at 
the various intersections of the fractures, each node being associated with a matrix 
volume, and determination of the flows between each fracture mesh and the associated 
matrix volume in a pseudosteady state. The method can be applied in hydrocarbon 
production well testing for example. 
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